/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  dev
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      testCalcNewton;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Find the root of f(x) = x + B*x/sqrt(1 + sqr(x)) - A
// using Newton-Raphson iteration

A       1;
B       1;

x       #calc
#{
    // Lookup the coefficients
    const scalar A = $A;
    const scalar B = $B;

    // Initial guess for x
    scalar x = 0;

    scalar x0 = x;

    do
    {
        // Store the previous iteration x for the convergence check
        x0 = x;

        // Temporary sub-function evaluations
        const scalar f1 = 1 + sqr(x);
        const scalar f2 = sqrt(f1);

        // Evaluate the function
        const scalar f = x + B*x/f2 - A;

        // Evaluate the derivative
        const scalar df = 1 + B/(f1*f2);

        // Update x
        x = x0 - f/df;

        // Test for convergence
    } while (mag(x - x0) > 1e-6);

    os << x;
#};


// Write a x-y data for plotting with
// gnuplot -e 'plot "fx.xy" us 1:2 title "f(x)" with lines; pause mouse close'
x0      -2;     // Start x
xn      2;      // End x
n       101;    // Number of x increments

y #codeStream
{
    code
    #{
        // Lookup the coefficients
        const scalar A = $A;
        const scalar B = $B;

        const scalar x0 = $x0;
        const scalar xn = $xn;
        const int n = $n;

        const scalar deltax = (xn - x0)/(n - 1);

        Field<scalar> x(n);
        Field<scalar> fx(n);

        forAll(x, i)
        {
            x[i] = x0 + i*deltax;
            fx[i] = x[i] + B*x[i]/sqrt(1 + sqr(x[i])) - A;
        }

        // Write the x-y data
        setWriter::New("raw")->write
        (
            ".",
            "fx",
            coordSet(true,"x", x),
            "f(x)",
            fx
        );
    #};

    codeInclude
    #{
        #include "setWriter.H"
    #};

    codeOptions
    #{
        -I$(LIB_SRC)/sampling/lnInclude
    #};

    codeLibs
    #{
        -lsampling
    #};
};

#remove y

// ************************************************************************* //
